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ABSTRACT 



In this paper we explore numerically the evolution of a warped accretion disc. 
While previous analyses have concentrated on the case where the disc is thick enough 
that the warp propagates as a wave, we focus here on the opposite regime of a thin 
disc, where the warp evolves diffusively. By comparing the numerical results to a simple 
diffusion model, we are able to determine the diffusion coefficient of the warp, a2, as a 
function of the relevant disc parameters, such as its thickness and especially its viscos- 
ity. We find that while in general the disc behaviour is well reproduced by the diffusion 
model and for relatively large viscosities the warp diffusion is well described by the 
linear theory (in particular confirming that the warp diffusion coefficient is inversely 
proportional to viscosity), significant non- linear effects are present as the viscosity 
becomes smaller, but still dominates over wave-propagation effects. In particular, we 
find that the inverse dependence of the diffusion coefficient on viscosity breaks down 
at low viscosities, so that a2 never becomes larger than a saturation value Ofmax of 
order unity. This can have major consequences in the evolution of systems where a 
warped disc is present. In particular, it affects the location of the warp radius in the 
Bardeen-Petterson effect and therefore the spin up (or spin down) of supermassive 
black holes in the nuclei of galaxies. Additionally, we also find that while the rate of 
warp diffusion does not depend significantly on the detailed viscosity formulation, the 
rate of internal precession generated by the warp is strongly affected by it. Such effects 
should be considered with care when modeling the evolution of warped discs. This em- 
phasises the need to test the above results using different numerical schemes, and with 
higher resolution, in order to investigate the degree to which numerical simulations 
are able to provide accurate modeling of the complex fluid dynamics of warped discs. 
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1 INTRODUCTION 

Many (perhaps most) astrophysically relevant accretion 
discs are expected to be twisted or warped. The torques 
responsible for such warping can be very different, in- 
cluding t idal effects of a c ompanion star in binary 
systems (jLarwood et al.l Il996l ). dynamical effects dur- 
ing the formation of protostellar discs, general relativis- 
tic Lense-Thirring precessi on ar o und a spinning black 
hole llBar dccn & Pcttcrsoj Il975l : IScheuer fc Feilerl 1 19961 : 
iLodato fc Pringld |2006|), and self-i nduced warping caused 
by radiation pressure (|Pringlel|l996l '). 

Observationally, warps are found both in discs which are 
expected to be relatively thick, such as the case of the hyper- 
accreting X- ray binary SS433 (|Begel man et al.ll2006D or, on 
the muc h less energetic side, the disc a round the young star 
KH 15D IChiang fc Murrav-Ciavll2004l '). In some other cases. 



warps are also fou nd in very thin discs, in the nucleus of the 
AGN NGC 4258 l|Herrnstein et al.lll996l : jPapaloizou et al. 
1998) or in the X-ray binary Her X-1 (jWiiers &: Pringle 

Ii99ay 

The secular evolution of such warped discs deeply af- 
fects the structure of the system. For example, the process 
of alignment between the spin of a black hole and the an- 
gular momentum of the disc due to the Bardeen-Petterson 
effect strongly dep ends on th e effectiveness of warp diffu- 
sion in thin discs (|Lodato fc Pringle ,2006, ). This can have 
major consquences on the spin up or spin down of the black 
hole during an accretion event, having implications also on 
the formation and ear ly growth of supermassive black holes 
iKing fc Pringlell2006l ). 

While the linear (and mildly non-linear) theory 
of warp prop agation has been in v estiga t ed exte n sively 
in the past (|Papaloizou fc Pringld 1 19831 : IPringld Il992i : 
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iPapaloizou fc LinI Il995l : lOgilviel 1 19991 . 12OO0I ) and has been 
tested numerically in the case where t he disc is thick enough 
that warp propagation is wave-like ()Nelson fc Papaloizoul 
[l999. . ■200Q '). there is as yet no numerical confirmation of the 
diffusive warp propagation expected in the thin disc case, 
nor a full description of the non-linear effects that might 
arise when the amplitude of the warp becomes large com- 
pared to viscosity. 

We consider here the propagation of warps in thin Kep- 
lerian accretion discs in the case when the disc is sufficiently 
viscous that the warp propagates in a diffusive manner. 

We start by defining the main disc properties. We con- 
sider a thin disc, rotating with angular velocity 0,{R), with 
surface density E(_R) and angular momentum per unit area 
L(i?). Here R should be interpreted as a 'spherical' coor- 
dinate. The local direction of L can be oriented anyhow in 
space, and the unit vector 1(_R) = L(7?)/L(i?) defines its di- 
rection. If the disc is rotating around a central point mass 
M, then its rotation is Keplerian, with Q, = 'J GM/R^ and 
L{R) = T,{R)VGMR. 

The disc is warped whenever the direction identified by 
1 changes with radius. We can characterize the warp ampli- 
tude with the dimensionless parameter tp, where 



tP = R 



dliR) 



dR 



(1) 



The disc thickness is H — Cs/f2, where Cs is the sound 
speed, and is the scale over which density and pressure 
change in the local z direction. The disc aspect ratio is H/R, 
and we shall assume that H/R <^ 1. 

IPapaloizou fc Pringld (|l983l ) and IPrindd (|l992[ l have 
introduced a simple equation to describe the evolution of 
a warped disc, in the case where the warp propagates di f- 
fusively in the disc. As shown bv IPapaloizou fc LinI (|l995h . 
diffusive behaviour occurs for isotropic viscosity if the disc 
aspect ratio H/R < a, where a is the viscosity parameter 
commonly used as a measure of the {R, (f)) (viscous) stresses 
in the disc (Shakura fc Sunyaev, 1973). In this case the warp 
prop agation can be approximately described by the equation 
l|Pringle„199a) 
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In this equation, the terms proportional to U\ describe the 
standard viscous evolution of a thin and flat disc, where ui 
is the usual viscosity acting on the [R, 0) shear. If we use 
the Q-parametrization, we have 



vi — aCaH = aD,H 



(3) 



The terms proportional to 1/2 in Equation ([2} arise 
whenever the disc is warped and |91/i9ii| 7^ 0. According 
to Equation ([2]) the warp diffuses with a diffusion coefficient 
1^2 ■ In view of this we shall define a second parameter 0:2 
defined similarly so that 



It is clear that the nature of the evolution of a warped ac- 
cretion disc is determined mainly by the relative values of a 
and 0:2. 

We stress that Equation ([2]) was not derived from any 
hydrodynamical (or magneto-hydrodynamical) equations. 
Rather the equation results simply from the conservation 
equations (mass and angular momentum), coupled with the 
assumption that the warp propagates in a diffusive man- 
ner. Thus there is no linear approximation required in the 
derivation and the equation is, in principle, applicable to 
warps of any amplitude, with an appropriate choice of 1^2- 
Thus, equation @ merely states that the warp evolution is 
diffusive, with a diffusion coefficient 1/2 to be determined. 

A number of authors have addressed the problem of 
warp diffusion. IPapaloizou fc Pringld l|l983l ) considered the 
internal hydrodynamics of the disc, assuming an internal 
isotropic viscosity for a warp small enough to be treated as 
a linear perturbation (i.e. <^ H/R), and for viscosity such 
that H/R < Q ^ 1. In this case the internal hydrodynamic 
flow is a laminar one, and is relatively straightforward to 
analyse (See Section 4). They found that to first order Equa- 
tion ((2]) gives the correct behaviour, along with the rather 
surprising result that 

(5) 



Q-2 



_ _ 

vi a 2o? 
This implies that 
1 



0.2 ~ 



2a 



(6) 



and therefore that the warp diffusion coefficient is inversely 
proportional to the size of the viscosity. We discuss the 
physics behind this in Section 3] To next order in a there 
are additional terms which imply precession of the disc, and 
which are not included in the simple formulation of Equa- 
tion 

lOgilvid (|l999l ) (see also lOgilvid[2000l : lQgilvie fc Dubud 
I2OO1F extends these approximate analytic results by use of 
an asymptotic expansion in terms of the small quantity 
H/R, but retaining the assumption of an isotropic (Navier- 
Stokes) viscosity. By this means he is able to take account 
of larger values of a and ip and to this approximation the 
internal dynamics still in the form of well-ordered, laminar 
flows. Ogilvio (1999) finds 



(7) 



1^2 



• a2CsH = a20,H 



(4) 



U2 _ 1 4(1 + 7q^) 
^ 2a2 4 -f a2 ' 

when ip <^ 1 and also gives higher order corrections with 
terms in i/;'^. He also finds additional precessional terms, 
which are generally of higher order than the terms already 
included (We discuss this further in Section [3.41) . 

None of these studies takes into account the fact that 
viscosity in real accretion discs is most likely not a sim- 
ple isotropic Navi er-Stokcs viscosity. It is generally believed 
(see for example iBalbus fc HawlevI Il998l ) that the mecha- 
nism responsible for angular momentum transport in discs 
(corresponding to i^i) is MHD turbulence, driven by the 
magneto-rotational in stability. An att e mpt at estimating 
U2 has been made by i Torkelsson et al.] (l2000f ) in a numer- 
ical shearing box simulation. It has become evident re- 
cently, however, that the results of such simulations are 
somewhat problematic ( King: et al]|2007l : iPessah et al ] |2007l : 
iFromang fc Papaloizoul 20071 ) . In addition, as discussed by 
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Pringle (1992), there is no particular reason to assume that 
the mechanism which damps the warp (corresponding to 
V2) is the same. In addition, for large (astrophysically sig- 
nificant) warps the internal motions can become sufficiently 
large, with relative supersonic velocities, that a number of 
authors have argued tha t they are hkely to be unstable 
IColeman fc Dopital [l993 : [ Gammie et al.ll2000l ) and so re- 
sult in enhanced dissipation. 

In the opposite regime of thick discs, or of lower viscos- 
ity, the w arp evolves to first order as a non-dispersive bend- 
ing wave (|Papaloizou fc Linlll995l ). The equations of motion 
for a wave in the case where the disc is Keplerian and ne arly 
inviscid are (|Lubow fc Ogilviell2000l : iLubow et aLll2002l ) 



(8) 
(9) 



— + anG = ^R n- — , 

where G is the disc inte rnal torque. The co rresponding dis- 
persion relation is then ([Nelson fc PapaloizouM,1999, ) 



iu{Lu — iaQ) — 



(10) 



where uj is the wave frequency, k the wavenumber and Cs the 
sound speed. The solution for uj is 



Lu=^{ian±[cik^ -a^n'']'-^^}. (11) 

Thus the above relation shows that in the absence of vis- 
cosity (a = 0) a bending wave propagates at half the sound 
speed. When a ^ Hk, i.e. when a ^ 2-nH/\, where A is the 
wavelength of the perturbation, we see that u is purely imag- 
inary, and wave-like propagation no longer occurs. Propaga- 
tion becomes purely diffusive in the limit \uj\ <^ aQ, in which 
case we obtain 



UJ 



f c,k y 

\2an J 



« 1, 



(12) 



Thus pure diffusion occurs when a > ^Hk = ttH/X. The 
boundary between the purely diffusive and purely wave-like 
propagation regimes is not a sharp one. For convenience 
we shall take the transition between diffusive and wave-like 
regimes to occur at the value a = Oc when the real and 
imaginary parts of ui are equal ljr = uji. Thus the critical 
value of a is given by 



V2 



Hk. 



(13) 



In this paper we investigate the hydrodynamics of dif- 
fusive warp propagation (i.e. a > Oc) using numerical sim- 
ulations. In particular, we use Smoothed Particle Hydrody- 
namics (SPH). Such a technique has been used in some of 
the mos t accurate investigations of warp propagatio n to date 
l|Nelson fc Papaloizoulll999i . l2000l : lLarwood et al.|[l99^ ) and 
is particularly suited to the treatment of complex and vari- 
able geometries, like that of a warped disc. However, as we 
will further discuss below, some numerical effects, for exam- 
ple related to the implementation of viscosity, might play 
a role and we stress the importance of conducting similar 
analyses with different codes, especially in view of simulat- 
ing a warped disc that evolves under the action of MHD in- 
stabilities, which are not taken into account here. The early 



< 
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Figure 1. Initial profile of the a;-component of the unit vector 
1 in our SPH simulations. Note that the y-component is initially 
zero, so that the initial warp has no twist. 



investigations bv lNelson fc Papaloizoul (|l999l . l2000 l'l have fo- 
cussed on the thick disc regime, where the warp propagates 
as a bending wave, but have already noticed some dissipa- 
tional effects. In this respect, the present work can be seen as 
complementary to these previous ones, since we here focus 
on the thin disc regime, where the character of warp propa- 
gation is expected to differ substantially from the thick disc 
case. 

The paper is organized as follows. In Section 2 we 
present the numerical setup, and we present the results in 
Section 3. In Section 4 we present a simplified discussion of 
the physics of diffusive warp propagation, and use this anal- 
ysis to interpret our results. We discuss the implications and 
limitations in Section 5. 



2 NUMERICAL SETUP 

We have performed a set of three-dimensional SPH simula- 
tions of warped discs. SPH is a Lagrangian hydrodynamic 
scheme (ILucvII 19771 : iGingold fc Monaghan|[l977l : lBenj|l990l : 



lMonaghanlll992l ). in which the gas disc is modeled with A'^ 



particles. For most of our simulations we have used 1 or 2 
million particles (these are some of the highest resolution 
SPH disc simulations to date). 

In the following we use code units in which the gravi- 
tational constant G = 1, the central point mass M = 1, so 
that at a radius i? = 1 in code units, the dynamical time is 
Q-^ = 1. 

We set up our warped disc by placing the gas particles 
in Keplerian orbits around a central poi nt mass M = 1 (i n 
code units), modeled as a sink particle l|Bate et al.lll995l ). 
with accretion radius R = 0.5 in code units. We distribute 
them in such a way that the disc has a prescribed initial 
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Table 1. Summary of the main physical properties of the simulations carried out. Here HjR refers to the value calculated at i? = i?o. 
See text (section [Sj for the definition of the parameter / and for details of the calibration of the disc viscosity parameter a. Note that 
OL2 = //2a. 




Figure 2. Azimuthally averaged vertical density profile for simulations SO to S9, at JJ = 1, where the disc thickness H = 0.02 (left) and 
R = 3, where H = 0.045 (right). The solid line shows the theoretical Gaussian profile expected in vertical hydrostatic balance, while the 
dashed and dotted lines (almost exactly coincident) show the density evaluated from the SPH code at t = 10 and t = 30 code units, 
respectively. Deviations from the expected profile only occur at a height z > 2H, where the density has decreased by a factor 0.2. 



surface density profile, as described below. We assign each 
particle a radius dependent sound speed Cs, so as to attain 
the desired temperature profile (see below). The particles 
are then distributed in z so as to attain a Gaussian density 
profile in the vertical direction, with thickness H — Cs/Q. 
We then tilt the orbit of each particle in such a way that 
the components of the unit vector 1 are given by 



( 

A 
2" 



1 + sin TT 



R2 -Ri 



for R< Ri 

for Ri <R<R2 (14) 
for R> R2 



where _Ri = 3.5 and R2 = 6.5 in code units, and Ro 



(15) 
(16) 



Warp diffusion in accretion discs 5 



(Ri + R2)/2. The initial shape of the warp is plotted in Fig. 
[1] The warp amplitude tp is then 



dl 



dR 



Rdh 
I, dR 



(17) 



the maximum of which is attained at R m Rg and is given 
by 

■kRoA 



2.62A 



2{Ri - Ri) 

The disc extends from Rin ~ 0.5 to _Rout 
surface density profile, E, given by 



(18) 
10, with a 

(19) 



Note that the density normalization So does not play a role 
in our simulations, since we do not include the disc self- 
gravity in the computation. The parameter p is set to p = 
3/2. The simulations we perform are locally isothermal, with 
a sound speed Cs, given by 




Ca,0-R 



(20) 



where 5 = 3/4 and where the normalization Cs,o determines 
the disc thickness. 

In order to model the viscous evoluti on of the disc, w e 
use the SPH artificial viscosity formalism (iMonaghanll 19921). 
It ha s been shown IjArtvmowicz fc Lubow 1994 : Murray! 
1 19961 ) that artificial viscosity in SPH can mimic the be- 
haviour of an Q-viscosity, with 



a OC QSPH' 



H 



(21) 



where qsph is the artificial viscosity coefficient, (h) is the 
average smoothing length and H is the disc thickness. Our 
setup has been chosen in such a way that the disc thick- 
ness is uniformly resolved at different radii. Indeed, the disc 
thickness H scales with radius as 



H 



OC R 



,3/2-g 



(22) 



and we expect the average smoothing length at one radius 
to scale as 

(h) OC p-^/^ OC f ii ] OC 



(23) 



Thus, by choosing p = 3/2 and q — 3/4, we have that 
H OC (h) OC J?^^''. This then produces a disc with a con- 
stant a, the magnitude of which can be adjusted by varying 
the SPH artificial viscosity parameter asPH- Additionally, 
since vi — aQH^, and H oc R^^^, we also expect that the 
simulations be characterized by a constant viscosity As 
for the quadratic term in the standard formulation of SPH 
artificial viscosity (the '/3'-term), we have only included it 
in simulations for which the 'viscous switch' is on (see be- 
low), in which cases we have adopted the common practice 
of setting /3sph = 2aspH. 

We have performed several simulations, varying the disc 
aspect ratio H/R, the viscosity a (through varying the SPH 
artificial viscosity asPH, see Section [3] below for details on 
how we calibrate the viscosity) and the peak warp amplitude 

^ We thank Gordon Ogilvie for suggesting this setup to us. 



?/). The main parameters of our simulations are summarized 
in Table [1] The ninth column of Table [T] shows the (radius- 
independent) ratio of the azimuthal and vertical average of 
the smoothing length h and the disc thickness H, which also 
demonstrates that we resolve the vertical structure moder- 
ately well in the simulations. Note that this column shows 
a vertical average of the smoothing length h, including the 
contribution of particles at high z which have a smaller den- 
sity and hence a larger h. The actual value of h in the bulk 
of the disc can be significantly smaller than the average. As 
a further check of the resolution achieved in our simulations, 
we plot in Fig. [5] the azimuthally averaged vertical density 
profile for simulations SO to S9 (see Table [l|, at two differ- 
ent radii, R — 1 (left panel) and R — 3 (right panel). The 
solid line shows the theoretical profile derived applying the 
standard hydrostatic balance condition, while the dashed 
and dotted lines (almost exactly coincident) show the value 
of the density as computed from the SPH scheme at two 
times, t = 10 (dashed) and t = 30 (dotted) code units. 
We can clearly see that the SPH estimates follow nicely the 
predicted profile, deviating from it only at vertical height 
z > 2H, where the density is already a factor 0.2 smaller 
than the peak. 

As can be seen from Table [T] we have considered three 
different values for the disc thickness. Most of our simula- 
tions have H/R « 0.0133, a few have H/R « 0.0334 and 
only one had a relatively large thickness of H/R ~ 0.0668. 
The critical value of a below which wave propagation is 
expected to occur is then Oc ~ 0.05, 0.12, and 0.25, re- 
spectively. It can be then easily seen that the thickest case 
we consider falls in the wave propagation regime, while the 
thinnest case (the majority of simulations) falls in the dif- 
fusion regime. The intermediate thickness is marginal, and 
wave propagation with significant dissipation is expected. 

2.1 On modeling disc viscosity in SPH 

In this paper, we have used the internal SPH viscosity in or- 
der to simulate a (Navier-Stokes) viscosity in the disc. This is 
clearly an idealization, since there is no a priori expectation 
that real disc viscosity in any way resembles Navier-Stokes. 
In practice accretion discs transport are likely to be be dom- 
inated by other forms of stress (for example, due to the 
magneto- rotational instability). But at least this provides 
some insight on the evolution of the warp in this specific 
case, and provides some comparison with previous analytic 
work. 

However, care should be taken when using the SPH ar- 
tificial viscosity to model a Navier-Stokes viscosity. Stan- 
dard implementations of SPH adopt the so-called 'viscous 
switch', so that viscous forces are only active when acting 
on approaching particles, and vanish for receding particles. 
Clearly, if we want to simulate something approximating a 
Na vier-Stokes viscosity we should tur n the 'viscous s witch' 
off. lArtvmowicz fc Lubowl (|l994l ') and iMurravl (|l996l ') have 
already remarked on this. 

Because it is unclear how the real, physical disc vis- 
cosity behaves, we have run simulations with and without 
the 'viscous switch'. The column 'VS' is Table [T] indicates 
the simulations that had the switch on or off (note that 
obviously in order to reproduce the same physical viscos- 
ity coefficient a, a much smaller value of aspu - roughly 
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Figure 3. Evolution of simulation Sll. The solid lines refer to the azimuthal and vertical average of the SPH simulation, while the 
dashed lines show the evolution of the corresponding initial conditions obtained by applying Eqs. ^ and ^ for the warp, while the 
surface density is evolved using a viscosity parameter a = 0.08. The left panel shows the evolution of the surface density S, while the 
right panel shows the evolution of Ix- The different lines refer to t = and to t = 155 (in units of the dynamical time at /? = 1). 



half - is needed when the viscous switch is oflt, as shown 
in Table [T]). This enables us to directly check the impact 
of using the viscous switch or not in the evolution of our 
warped discs. As we discuss below (Section 3.3), the effects 
on the evaluation of the warp diffusion coefficient (i.e. on the 
diffusion/damping of the warp) are minor, if not negligible. 
There is however, a strong sensitivity in the rate, and indeed 
direction, of internal precession induced by the warp. 



3 RESULTS 

In this Section we describe the main results of our simu- 
lations. Our general procedure for the analysis is relatively 
simple. We first divide the radial range of the simulation in 
a series of annuli at different radii R. From the output of the 
SPH simulation we then compute the surface density E(i?) 
and the azimuthal and vertical average of the three com- 
ponents of the unit vector 1 as functions of time. We then 
compare the evolution of these quantities with the evolution 
predicted either by Eq. ([2]) (for the thinnest cases), or by 
Eqs. (O and ^ (for the thick cases). 

3.1 Wave-like propagation in the thick disc case 

As mentioned above, we expect the warp to propagate as a 
bending wave in the thickest of our simulations, and possi- 
bly also the intermediate thickness ones. This is indeed con- 
firmed by our simulations. As an example, we show in Fig. [3] 
the evolution of simulation Sll. The left hand panel shows 
the evolution of the surface density S (solid lines), compared 
to the evolution predicted by a simple viscous surface den- 
sity evolution with a — 0.08. The right hand panel shows 



the evolution of the x-component of 1 (solid lines) compared 
to the evolution predicted from Eqs. ^ and @. Note that 
in order to obtain a good fit to the simulations it is essen- 
tial also to include the dissipative term in Eq. ((9}, i.e. the 
second term on the left-hand side of the equation, since a 
pure wave evolution does not fit the data. We also note that 
the same evolution could also in principle be fitted (at least 
over the limited time span covered by the simulation) with 
a diffusive warp propagation model, but where the warp dif- 
fusion coefficient is taken to be much smaller than the value 
predicted from the linear theory, i.e. Q2 « 2.5 ^ l/2a. 



3.2 Diffusive propagation in the thin disc case 

For our thinnest disc case we expect the warp to propa- 
gate diffusively, according to Eq. ((2}. We have then solved 
Eq. ^ numer i cally using t he techniques desc r ibed in de- 
tail in lPringlel (| 19921 ) and in iLodato fc Pringld (|2006D . The 
model described in equation ((2} contains two free parame- 
ters, namely, the two viscosities and 1/2, and we adjust 
these in order to match the evolution of the simulation. In 
particular, the surface density E is most sensitive to the 
value of the viscosity z^i, while the evolution of the compo- 
nent Ix of the unit vector 1 is most sensitive to the value 
of V2- We therefore use these two quantities separately to 
determine the two paramet ers. To be su re, the evolution of 
E does also depend on U2 (|Pringl6lll99"3 ). However, this ef- 
fect scales quadratically with the warp amplitude and can 
be generally neglected in most of our simulations (with some 
exceptions, as shown below). 

In practice, we determine the value of a and the value 
of a parameter /, which is a measure of the deviation of 1/2 
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Figure 4. Evolution of simulation S2. The solid lines refer to the azimuthal and vertical average of the SPH simulation, while the dashed 
lines show the evolution of the corresponding initial conditions obtained by applying Eq. l[2]l, with the following parameters: a = 0.18 
and / = 1. The left panel shows the evolution of the surface density S, while the right panel shows the evolution of l^- The different 
lines refer to t = and to t = 465 (in units of the dynamical time at i? = 1). 




Figure 5. Evolution of simulation S6. The solid lines refer to the azimuthal and vertical average of the SPH simulation, while the dashed 
lines show the evolution of the corresponding initial conditions obtained by applying Eq. l[2]l, with the following parameters: a = 0.07 
and / = 0.42. The left panel shows the evolution of the surface density S, while the right panel shows the evolution of l^. The different 
lines refer to t = and to t = 550 (in units of the dynamical time at = 1). 
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Figure 6. Evolution of simulation S6. The solid lines refer to 
the azimuthal and vertical average of the SPH simulation, while 
the dashed lines show the evolution of the corresponding initial 
conditions obtained by applying Eqs. ^ and JSJ, with a = 0.07. 
The different lines refer to t = and to t = 550 (in units of the 
dynamical time at R= 1). 

from the value expected from the linear theory. Thus / is 
defined by 



The fifth and sixth columns of Table [T] show the mea- 
sured values of a and / for all of our simulations. It can be 
seen that we are able to span a range of a that goes from 
0.05 to 0.28, that is roughly a factor 7. Only for the three 
largest values of a does the measured value of / agree with 
the linear estimate, while for most other cases the diffusion 
coefficient U2 is smaller than predicted by linear theory. 

Note that, with regard to the evaluation of the warp 
diffusion coefficient, the simulations which have an imple- 
mentation of viscosity that more closely resembles a Navier- 
Stokes viscosity (i.e. viscous switch off) do not show sig- 
nificant differences with respect to the corresponding ones 
with the viscous switch turned on. Simulations SO and SOb 
and SI and Sib, that have the same estimated value of a, 
also share the same warp diffusion coefficient. Simulations 
S3 and S4, which have very similar values of a, also have 
a similar diffusion coefficient. Note that this agreement en- 
compasses the whole range of viscosities that we probe in 
the thinnest disc configuration, extending from high values 
of Of, for which the warp diffusion rate agrees with the linear 
theory, down to lower values, where some non-linear effects 
start to play a role. 

Fig. IHshows one example of our analysis, corresponding 
to the case where H/R — 0.0133, a = 0.18 and i/imax ~ 0.026 
(simulation S2). The left panel shows the evolution of E, 
while the right one shows the evolution of Ix . The solid lines 
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Figure 7. Evolution of simulation S8. The solid line shows the 
results of the SPH simulations, while the dashed line shows the 
diffusive evolution, according to Eq. ((2J. The two snapshots refer 
to i = and to t = 870 in code units. The inflection in the surface 
density is c aused by the e ffect of 1^2 on the surface density, as 
discussed bv lPringlj 1 I1992I) . 

refer to the results of the SPH evolution at two different 
times, t — and t = 465 in code units (where the unit time 
is the dynamical time ai R — 1), while the dashed lines show 
the evolution of the simple model of Equation ([2} . As can be 
seen, the evolution of the two quantities is well reproduced 
by the model, therefore demonstrating numerically for the 
first time the validity of the diffusive model for warp propa- 
gation in this regime. In this particular case, not only does 
the warp evolve diffusively, but the value of the diffusion co- 
efficient agrees with the expectations from the linear theory, 
that is, we measure / = 1. 

Fig. [5] shows the evolution of simulation S6. This sim- 
ulation is also well fitted by a diffusive model, except that 
in this case, although the value of a = 0.07 is above the 
critical value so that diffusive propagation is expected (in 
this case Uc — 0.05), the required value of Q2 ~ 3 is sig- 
nificantly below the value expected from the linear theory 
l/2a ~ 7.14. To be sure that wave propagation effects do 
not play a significant role here, we have also tried to fit the 
evolution of this simulation with a wave propagation model, 
as done in the previous section. The results are shown in 
Fig. |6] It can be seen that in this case a wave propagation 
model (including the effects of viscosity) does not reproduce 
the results of the simulation. 

Most of our simulations have a low amplitude warp, 
with ^ ~ 0.026. However, we have also run some simula- 
tions with ~ 1. In the thin disc case, these simulations are 
S8 and S9, for which ^ « 1.3. Simulation S8 has a — 0.26, 
for which the corresponding low amplitude simulations fol- 
low the linear predictions for a2 (Eq. (|6])). Simulation S9 
instead has a = 0.1, for which the low amplitude simulation 
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is already in the saturated regime for 02. In both cases the 
diffusion coefficient is found to be smaller than the value 
predicted from Eq. ((6| and we measure / — 0.75 for S8 and 
/ = 0.6 for S9. The results for simulation S8 thus show that 
increasing the warp amplitude leads to a marginal reduction 
of the diffusion coefficient. On the other hand. In the case of 
simulation S9, with q = 0.1 for which a low amplitude warp 
is already enough to induce a reduction of 0:2 with respect 
to linear theory, a further increase of does not lead to a 
significant further reduction of 0:2. 

These large amplitude simulations are also interesting 
because they reveal the effects of the warp diffusion term on 
the evolution of the disc surface density. Fig. [7] shows the 
evolution of simulation S8 (solid line) along with the evolu- 
tion of the simple diffusive model of Eq. ([2| (dashed line). 
Note how the model reproduces well the inflection of the 
surface density close to the warp location. The inflection is 
due to the effect of 1^2 on the evolut ion of E, as mentioned 
earlier, and shown in IPringld (|1992| ). In contrast with the 
previous low amplitude simulations, in this case, where the 
warp amplitude is large, this effect docs play a role. Such 
effect has b een described clear l y in[Pringle (,19921 and more 
recently in iLodato fc Pringld (|20oi)7 Note that such dissi- 
pative effects cannot be reproduced by a purely wave-like 
warp propagation. 



3.3 The diffusion coefficients 

Each one of the high resolution simulations described in the 
previous section is quite demanding in terms of computing 
time, and we could only run a few such simulations. As a 
result, we could only obtain a partial coverage of the relevant 
parameter space. Even so, we can still draw some interesting 
conclusions from the limited available data. 

First we find that the evolution of the disc is well de- 
scribed by Equation ([2}, and the numerical results can be 
fit by just varying the parameters a and 02 . We stress that 
we do not perform an actual statistical fit of the viscosity 
coefficients, but simply choose them so as to match the evo- 
lution of the numerical simulations. We estimate that the 
typical uncertainty on these parameters is « 7% for a and 
Ri 10% for a2. 

Fig. |S] shows the values of the diffusion coefficient for 
the warp, as obtained from all of our thin disc simulations, 
for which propagation of warp occurs diffusively. The solid 
symbols refer to the case where the viscous switch is turned 
off (that should be closer to a Navier-Stokes viscosity) , while 
the open symbols refer to the case where the viscous switch 
is on. For the two simulations with the highest a, the point 
refers to both cases, as the results were identical. The cir- 
cles refer to the high-amplitude simulations. As an example, 
we also show for one point the typical uncertainties of the 
estimated parameters. 

It can be seen that the points follow the expected rela- 
tion 0:2 = 1/2q: at large a, but as a decreases, the value of 
a2 begins to deviate from the theoretical relation and ap- 
pears to saturate at a value of around Ofmax ~3 - 4. However, 
perhaps surprisingly, comparison of runs S5 and S9 shows 
that we do not observe any dependency of the saturation 
value on the warp amplitude. 



Figure 8. Results of the numerical simulations. The points indi- 
cate the values of the diffusion coefficient 02 as a function of the 
viscosity coefficient a. The solid symbols refer to simulations that 
do not use the 'viscous switch' and should thus have a viscosity 
more closely approximating Navier-Stokes. The open symbols do 
use the viscous switch. The squares refer to the small warp am- 
plitude case Vmax ~ 0.026, while the circles refer to the large 
amplitude case Vmax ~ 1.3. The error bars shown represent the 
typical uncertainties on the diffusion coefficients. The solid line 
shows the expected value of 02 from the linear theory. It is ev- 
ident that our simulations reproduce the expected results from 
linear theory for values of q > 0.16. Below this value we find that 
02 appears to saturate at a value around Omax 3 - 4. 



3.4 Precession 

Our simulations also display some small precessional effects. 
This is shown in Fig. |9] that shows the evolution of the y- 
component of the angular momentum of the disc, ly{R). This 
is initially set to be zero (see Eq. (|15p l but it then grows 
due to internal precessional torques. Such torques are not 
accounted for in the diffusion model bv [Pri nglc (19 92J, b ut 
they are present in the full linear theory fe.g. fOgilviel ljl999h '). 
In interpreting our results we have added such terms in our 
simple diffusion model, by adding a ter m on the righ t-hand 
side of equation ((2)), in the form of (see IOgilvielll999l ) 



dt 



(25) 



where we have introduced a third coefficient 1/-^ related to 
precessional effects (with a corresponding 03 = vs/flH^). 
The solid lines in Fig. [9] show the evolution of the SPH sim- 
ulations, while the dashed lines refer to the simple model 
(including precession). The two simulations shown in Fig. 
[9] are Sib (left panel) and SI (right panel) which have the 
same a and 0:2 but use a different implementation of viscos- 
ity, so that SI uses the viscous switch, while Sib does not. 
The warp evolution of these two simulations agrees with the 
predictions of linear theory for the evolution of the warp. 
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Figure 9. Evolution of the component ly of the unit vector 1. Initially ly = and it is the precessional torques, neglected in Equation 
(2), which cause ly to evolve. Left panel: simulation Sib (with the viscous switch turned off). Here the precession occurs in the direction 
predicted by the linear theory, but with a smaller precession rate. The solid line represents the results of the SPH simulation, while the 
dashed line represents the result of a simple diffusion model including a precessional term, with 03 ^ 0.17. The curves refer to t = 825 
(in units of the dynamical time at i? = 1). Right panel; same plot for simulation SI, which has the same a and 02 as Sib, but the viscous 
switch is on. Here precession occurs in the opposite direction and the best diffusion/precession model able to reproduce the results has 
03 Si —0.46. Here the curves refer to t = 840 (in units of the dynamical time at R = 1). 



They do not, however, agree with linear theory with regard 
to disc precession. Moreover, there is an obvious difference 
between the two cases in that precession occurs in different 
directions. In order to model the evolution, we therefore re- 
quire a negative Ofa term in the case of SI and a positive 
Q3 for Sib. Note that the linear analytic theory, using an 
isotropic Navier-Stokes viscosit y predicts a positive with 
magnitude g iven by Q3 =3/4 jPapaloizou fc Pringle|[l983 : 
IOgilvielll999l ). The resulting values of Q3 for the simulations 
where we measure precession are shown in Table [1] None of 
our simulations agrees with the analytic results. 

Note that in all cases 03 is much smaller than 02, in- 
dicating that precession occurs on a much longer timescale 
than the timescale for diffusion and/or damping of the warp. 
This implies that, since internal precession is induced by the 
warp, it is unlikely to play any role in the disc dynamics, be- 
cause by the time it starts to be important, the warp has 
already been diffused out, at least in the case where there is 
no strong external torque to enforce a given warped shape. 



4 DISCUSSION 

In this Section we discuss the outcome of our numerical sim- 
ulations and compare them with the predictions of the ap- 
proximate analytic theories. Our general finding is that for 
values of a > 0.16, the value of 02 follows the simple theo- 
retical prediction that a2 ~ l/2a. However, for lower values 
of a our simulations indicate that the value of 02 seems to 
saturate at a upper limit of 

Q^2 — Q^max 4. 



4.1 The theory of the simple Q2 a relation 

In order to aid theoretical understanding of the processes 
involved, we first introduce a simple, physical model to de- 
scribe the diffusive evolution of a warp in thin, viscous discs. 
This model essentially reproduces the more rigorous analysis 
of Papaloizou & FriuElc (198311 . describing the main physi- 
cal processes in a simple and intuitive way. 

The evolution of the warp is controlled essentially by 
hydrodynamical processes taking place inside the disc. The 
main point is that a warp, coupled with the vertical strati- 
fication (with typical scale-height H) in a thin disc, induces 
a horizontal pressure gradient which depends on vertical 
height z (see Fig. llOp . Thus, as a fluid element orbits in the 
disc with a frequency f2, it experiences a horizontal (radial) 
pressure force, which oscillates at the same frequency Q,. The 
crucial point is that for a Keplerian disc this forcing term is 
resonant with the epicyclic motion, since the epicyclic fre- 
quency K = Q. for a Keplerian disc. Hence the presence of 
the warp excites strong horizontal epicyclic motions within 
the disc. The magnitude of the induced pressure gradient is 
given roughly by 



dp 



V" 



pi/) 

H • 



(26) 



dp 

Thus the corresponding force term in the horizontal momen- 
tum equation is 



1 ap 

pdR 



clip 
H 



(27) 



where we have also used p ^ c^p. 
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Figure 10. Schematic view of a warped disc. The shaded areas 
indicate regions of higher pressure. The arrows indicate the di- 
rection of the horizontal pressure gradients induced by the warp. 
As a fluid element orbits around the centre, it feels an oscillating 
radial pressure gradient, whose amplitude is a linear function of 
the height z. 



As is shown in more detail in Papaloizou & Lin (1995) 
the forcing term varies linearly with z, and therefore excites 
a response for example in the radial velocity field which takes 
the approximate form 



Vr 



: Vj_ 



(J)cos0. 



(28) 



Here is the azimuthal angle in the disc, and vj_ denotes the 
induced horizontal velocity in the radial direction at the disc 
'surfaces' z = ±H. If we let x{t) be the radial Lagrangian 
displacement of a disc particle at the disc surface relative to 
an unperturbed disc particle on a circular orbit, then this 
implies that 



x{t) 



n 



sin fit. 



(29) 



Then the disc particle at position x{t) is subject to three 
forces: a restoring force which gives rise to the epicyclic fre- 
quency K = n, a damping force due to some form of vis- 
cous dissipation, and a pressure force caused by the warp 
with amplitude given by Equation (|27p . Thus x{t) obeys an 
equation of the form 



(30) 



x + - + n'^x = Rn'^ (^■^^ ™^ 

Here the amplitude of the forcing term on the r.h.s. is taken 
from Equation (|27|) . The third term on the l.h.s. is the 
epicyclic restoring force. The second term on the left-hand 
is a damping term, where we have introduced a damping 
time-scale r, that describes the viscous damping of the hor- 
izontal shear motion (described in Equation I28|l induced by 
the warp. If the horizontal shear motion is damped by a 
kinematic viscosity then we have simply that 



(31) 



(Note that the typical lengthscale of the horizontal shear 
is H, rather than R, hence the difference between the 
above formula and the corresponding one associated with 
the viscous timescale in a differentially rotating disc, fvisc ~ 
RVi^i). 

If we make the further assumption (which, as we have 



stressed, is not necessarily true for a disc) that the viscosity 
is isotropic, so that Vz ~ i^i, we obtain 

- _L 



(32) 



In line with our expectations (Equation [29} we look for 
a solution to Equation ([30} in the form 



x{t) — asinQt, 



(33) 



where now a is the amplitude of the horizontal motions at 
the disc surfaces z = ±H. Substituting ([33} into (|30|) . we 
find a solution for the amplitude of the horizontal motion a 

a = ijHinr), (34) 

which corresponds to a horizontal velocity of 

v± — aQ — i/;Cs(r2r). (35) 

The amplitude of the horizontal shear motion, dva/dz, is 
then given by 



H 



(36) 



If we add the assumption that the horizontal shear is 
damped by an isotropic viscosity with magnitude measured 
by the usual a, this, using Equation (I32p . then implies that 



a ~ —H, 

a 



and the horizontal velocity is then given by0 



Vl^ = — Cs 

a 



(37) 



(38) 



With this assumption the amplitude of the vertical 
shear, S = dvji/dz, is given by 

S^^ = U (39) 

In order to relate the above description with the model 
described by Equation ([2}, we note that in that equation 
the coefficient z/2 acts formally on the radial derivative of 
the vertical component of the velocity (since it dissipates 
the warp). By the argument described above, the dissipation 
actually comes about because of damping of the resonantly 
induced, horizontal, shear motion, that communicates the 
warp between different radii. The definition of 1/2 inherently 
implies that the rate at which energy is dissipated in the two 
processes must therefore be the same. Thus we must have 



(40) 



where angled brackets imply suitable vertical and azimuthal 
averages. 

Now, in a warped disc Vz = i/'-Rf2 and du^/di? = tp^l, 
which, in equation (|40|) gives 

(41) 



U2 = uzi^iry = nn'^inr) 



where, in the last equality, we have used the definition of Vz, 
Equation (|31|l . This is equivalent to writing 



^ As an aside, we emphasise that this estimate is not valid for 
small values of a, because for a < ac the propagation becomes 



wavelike, and in that case v 1 



(l/'/Oc)Cs 
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Q2 = ^T. (42) 

If the damping time-scale r is indeed given by isotropic 
viscous dissipation, then, as above, we have tQ. = 1/a and 
we finally recover the desired scaling of V2 with a in the 
linear case (cf. equation ([5j) 



It is worth noting here that this result is somewhat 
counter-intuitive, in that it implies that the damping of 
the warp occurs more rapidly when the viscosity is smaller! 
This comes about because a small value of isotropic viscosity 
(small a) permits a large resonant velocity {vr, oc 1/a). The 
dissipation rate depends on vv\ oc l/a, and so increases as 
Q decreases. It is evident that the analysis must break down 
at some stage as a ^ 0. 

4.2 Comparision with the simulations 

Most of the results shown in Figure |S] come from a set of 
simulations all with H/R = 0.01333 and i/imax = 0.026 but 
with varying values of a. For these simulations we find that 
relationship between a and 02 predicted by simple linear 
theory holds for values of a above around a ~ 0.16 where 
02 ~ 3. Below this critical value of a, 02 remains approxi- 
mately constant. 

We now need to ask what mechanism(s) might give rise 
to such an effect, and in addition whether they correspond 
to some valid physics, or are merely numerical artifacts. 

We have noted that from simple physical considerations 
we expect that 02 ~ Sir, where r is the timescale for damp- 
ing the warp-induced horizontal shearing motions. Thus an 
upper limit to a2 can come about through there being an 
upper limit to the damping timescale r; that is, we need 
to identify a mechanism which prevents the damping tak- 
ing place too slowly. For isotropic viscosity, the damping 
timescale is such that Q.T — 1/a, so that the lower the vis- 
cosity, the longer the timescale. Thus, paradoxically, in order 
to set an upper limit on 02, i.e. in order to stop the warp be- 
ing dissipated too fast, we need to find a mechanism which 
gives rise to increased dissipation. If there is such a process 
such that r < Kfl^^, where K is some constant, then we 
would find that a2 < K — Ofmax, in line with our findings. 

So what physical or numerical process(es) might give 
rise to enhanced dissipation for small values of a? 

We have noted that for isotropic viscosity the Mach 
number of the horizontal velocity induced by the warp is (to 
within factors of order unity) 

vjl/cs =M = S/Q, = ip/a. (44) 

The above equation shows that the horizontal resonant mo- 
tion becomes supersonic for a ^ ■;/). In our simulations 
breakdown of the simple relation occurs for a ~ 0.16 and 
tp ~ 0.026, and therefore when M « 0.16. Thus, breakdown 
occurs when the velocity difference across the disc thickness 
\v^{H)-vj_{-H)\ « 0.32c, 

The radial shear flow induced by the warp can be clearly 
noted in Fig. 1111 which displays the radial Mach number for 
simulation S6 at t = 195, at i? = 5 and at different azimuthal 
positions, that is cjf> = (upper left), <j) = 11/2 (upper right), 
(j) = TT (lower left) and <j!> = 37r/2 (lower right) (Here, we 



use a cylindrical coordinate system such that the z-axis lies 
along the local direction of the angular momentum vector 
1). At height z = ±H the amplitude of this radial oscillation 
is v±/cs ~ 0.1, much less than the value of 0.37 expected for 
this simulation based on Eq. (|44|l . 

We have also run two simulations (S8 and S9) with a 
much larger warp amplitude, that is ^max ~ 1.3. For sim- 
ulation S9, we see that a ~ 0.1, in which case we expect 
the Mach number of the radial resonant motion to be of the 
order of ~ 13, which is much larger than the Mach num- 
ber A4 « 0.16, above which we have observed the standard 
relation between a and a?2 to break down. This simulation 
should be compared with simulation S5, which has a very 
similar a and a very similar 02, but a much smaller tp and 
therefore a smaller predicted Mach number of TVf ~ 0.29 
(but still large enough to be in the saturated regime). We 
then see that increasing further the warp amplitude does 
not produce a signiflcant decrease in 02 . Thus although for a 
larger warp the amount of dissipation increases, the dissipa- 
tion timescale appears to remain unchanged at r ~ 3$!"^. In 
the case of simulation S8, instead, we have a ~ 0.26, so that 
the predicted Mach number would be « 5. The corre- 
sponding low amplitude simulation (SO), which has a ~ 0.28 
and predicted A4 « 0.1, lies in the non-saturated regime, 
and in this case increasing the warp amplitude does result 
in an increased dissipation and a reduction in the value of 
Q2, which for S8 lies marginally below the predicted relation. 

4.3 High Mach number and enhanced dissipation 

We have established that the results of our simulations can 
be understood if the internal disc flows, induced by the warp, 
result in higher dissipation than would be predicted by linear 
theory (assuming Navier-Stokes viscosity) when the Mach 
number of the flows exceeds a value of around M ~ 0.16, 
i.e. when the relative flow velocity across a full disc thick- 
ness exceeds around 0.32cs. The internal flows predicted by 
the analytic theory are in fact of two types (see, for exam- 
ple. Equations (102) and (103) of Ogilvie, 1999). The first 
is the sinusoidally oscillating, linear shear flow discussed in 
Section |4]abov43, which takes the form vr cc z and oscillates 
with frequency Q.. This flow has an amplitude which scales 
linearly with and in the linear theory is indeed propor- 
tional to xp/a. The second, whose magnitude depends more 
strongly on the size of the warp (that is, it scales with tp^), 
is a sinusoidally pulsating, homologous flow, which is alter- 
nately expanding and contracting with angular frequency 
2f2 and is of the form cc z. We have calculated the verti- 
cal gradients of these two types of flows by fitting a straight 
line to the velocity data, such as those shown in Fig. 1111 We 
plot in Fig. [T2]the results of this procedure as a function of 
azimuth for simulations SOb (left) and S9 (right). The typ- 
ical uncertainty for these data points ranges from 5 to 10 
per cent. In both plots the solid line refers to the gradient of 
Vr, while the long-dashed line refers to the gradient of v^. 
In the sob case, the short-dashed line is a sinusoidal fit to 
the data points for vr, with an amplitude equal to 0.063. 

^ This radial {vr) flow is part of a resonant epicyclic motion (see 
Section 4.1). It is therefore accompanied by an azimuthal flow of 
the form on z, with amplitude vr = 2v^. 
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Figure 11. Run S6. The radial Mach number induced by the warp as a function of height above the disc at i? = 5 and (/> = (upper 
left), </< = 7r/2 (upper right), <j) = it (lower left) and (p = 3tt/2 (lower right). Here the cylindrical coordinate system is such that the z-axis 
lies along the local direction of the angular momentum vector 1. The disc scaleheight here is -ff = 0.0665. The amplitude of this flow at 
z = ±H is 0.1. 



Let us first consider the case of simulation SOb. As dis- 
cussed above, this simulation agrees well with the linear the- 
ory and with its prediction for the warp diffusion coefficient. 
In agreement with linear theory, for this low amplitude warp, 
we find that the induced resonant flow is dominated by the 
radial component, the vertical one being negligible. The ra- 
dial flow is indeed oscillating with frequency Q as predicted 
by Ogilvie (1999). We have shown above (Eq. ^) that 
the amplitude of this oscillation should be (within factors 
of order unity) of the order of i/)/q ^ 0.093 for SOb. Based 
on the results shown in Fig. [12] we can estimate the exact 
proportionality factor to be 



dz 



0.65-0. 



(45) 



We can now check the predicted scaling with t/j/o as we 
decrease the value of a. To this end, we have repeated 
the above analysis for S4 and S6, and we have found in 
both cases a roughly sinusoidal oscillation of {dvii/dz)/Q,, 



with amplitude equal to 0.1 — 0.54i/)/a for S4 and to 
0.12 = 0.33'i/)/a for S6. Hence, we do not conflrm the pre- 
dicted scaling, the value of dva/dz for S6 being a half of its 
predicted value. This is further evidence for enhanced dissi- 
pation. Note that this enhanced dissipation is only found for 
small values of a and it appears to increase as we increase 
the value oi xp/ot. 

The presence of both the radial and the vertical flow 
can be seen clearly in the right panel of Fig. [12] and in Fig- 
ure 13 which corresponds to the run S9 which has a strong 
warp -i/) ~ 1. The structure of this pulsating flow is shown 
in Figure 14, where we show the Mach number in the z- 
direction as a function of z at various azimuths in the disc 
(as before, the cylindrical coordinate system used here is 
such that the z-axis lies along the local direction of 1). As 
can be seen, the flow is highly supersonic, with disc ele- 
ments approaching each other at velocities up to « lOcs . As 
shown by Ogilvie (1999) these flows are formal solutions to 
the analytic equations, even for large values of when the 
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Figure 12. Vertical gradient of the resonant flow induced by the warp for simulations SOb, for which i/i 0.026 (left) and S9 (right), for 
which -0 ^ 1.3. The linear theory predicts to first order in ip that Avu/dz ex. ■(/)sin(0 — if>o) and dv^/dz oc i/"^ sin2(</) — 0o)- The solid line 
shows dvji/dz, while the long-dashed line shows dv^/dz. The short-dashed line is a sinusoidal fit to dvn/dz for SOb. Here the cylindrical 
coordinate system is such that the z-axis lies along the local direction of the angular momentum vector 1. 



Mach numbers of the flows are expected to be large. Ogilvie 
(1999) does not, however, discuss the stability of such flows. 
Note that while both these flow appear to have the cor- 
rect azimuthal periodicity (with frequency SI for the radial 
component and 2Q for the vertical one), they are not sinu- 
soidal as additional Fourier components become important 
for large ^ (Ogilvie 1999, Appendix). Note, in particular, the 
large amplitude and the sharp jumps of the vertical compres- 
sion. Across a very narrow azimuthal range (and therefore 
over a timescale much shorter then dynamical) the velocity 
changes by a factor much larger than the sound speed. Given 
the short timescale, this change can be only due to pressure 
forces and artificial viscosity in this strongly compressional 
motion. Note also that the negative amplitude of the com- 
pressional motion is larger than the positive one, indicating 
that a large amount of kinetic energy is dissipated in the 
process. This also implies that the pulsating flow is not of 
the form v, (x z at all times. 



4-3.1 A numerical explanation? 

One possibility we need to address is whether the enhanced 
dissipation that we find in our simulations at high Mach 
number is a purely numerical effect. From the simulations 
presented here it is hard to be certain that it is not. But it 
is not easy to see what the effect might be. 

It is certainly true that the fiows are only moderately 
well resolved in our simulations: in order to treat discs which 
are thin enough that warp propagation occurs diffusively we 
were only able to consider simulations in which the aver- 
age number of SPH smoothing lengths across the thickness 
of the disc is around 3.3. But the explanation cannot sim- 
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Figure 13. Colour plot of the vertical density structure of simul- 
tion S9, at t = 185, close to the warp location at _R = 5. The 
arrows indicate the velocity structure. The vertically shearing ra- 
dial flow induced by the warp has, in this large warp amplitude 
case, a signiflcant component along the local vertical direction, 
determining a strong compression and rarefaction of the disc. 



ply be due to lack of resolution, since our simulations agree 
with the analytic predictions for for the larger values 
of a, in the regime where the induced velocities are low, 
and analytic solutions are expected to be a good approxi- 
mation. Disagreement only occurs once the Mach number 
of the induced shear flow increases to about 0.16 across a 
scaleheight H, or to about 0.096 across a smoothing length 
h. In addition we obtain similar results for a2 whether the 
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SPH viscous switch is turned on or off. Thus it seems to us 
unlikely that either lack of resolution or the nature of the 
SPH stress tensor is the cause of the disagreement. 

It is evident, however, that the disagreement between 
the approximate analytic theory and the numerical simu- 
lations in respect of precession rates does depend on the 
nature of the numerical viscosity. Our finding that none of 
our simulations confirms the analytic predictions, and that 
even for the same formal a we obtain opposite signs for the 
precession rates depending on whether or not the viscous 
switch is on or off, implies that the induced precession is 
very sensitive to the precise details of disc viscosity. This 
implies that if the viscous stress tensor within the disc is 
not precisely of the Navier-Stokes form, and as we have re- 
marked above it is very unlikely that it is, then the analytic 
predictions for the precession rates are unlikely to be rele- 
vant to real accretion discs. Luckily, as we have remarked 
above, such precession is in general only a minor effect. 

4-3.2 A physical explanation? 

When the Mach numbers of the internally induced flows ap- 
proach unity it seems physically plausible that such flows 
might become unstable. If such instabilities gave rise to en- 
hanced dissipation, then they could provide an explanation 
for our results. 

Instabilities in a steady shear flow typically grow at a 
rate 7 ~ S, the rate of shear (e.g. lOrazin fc Reidlll98ll) . In 
an oscillating shear flow one might expect instability to be 
suppressed unless the growth rate exceeds the rate of oscilla- 
tion i.e. unless S > Q. However, Gammie et al (2000) demon- 
strate that shear instabilities, caused by parametrically un- 
stable modes, can exist with 7 ~ 5 even when S < (and in 
particular they predict instability for S'/fi ~ > 30a). For 
the largest scale modes of size ~ H the damping timescale 
would be ~ H^/v ~ r ~ (aO,)'^ . Thus we might expect 
shear instabilities to be able to grow when 7r > 1, that is 
when -i/) > or when M > a. Using these arguments, for 
the set of simulations with tp — 0.026, we would then expect 
shear instabilities to occur when a < 0.16, which is in agree- 
ment with our findings. To discover whether this agreement 
is fortuitous or not will require further work. 

We have not demonstrated, however, that a shear flow 
becomes unstable simply because it is supersonic, and to 
this extent the arguments we put forward here are specula- 
tive. Indeed one could argue in this case that the particles 
on the top and bottom of the disc are essentially follow- 
ing free, elliptical trajectories around a central point mass. 
However, this is not strictly true because fluid elements sep- 
arated vertically by ~ _ff do need to interact, and to receive 
a vertical velocity nudge of order ~ Cs twice per orbit, i.e. 
on a timescale ~ (1/2)57"^, in order to remain a distance 
~ H apart. Once the horizontal fluid motions become su- 
personic, i.e. once S ~ il, the fluid elements separated ver- 
tically by ~ _ff can no longer maintain pressure communica- 
tion. Then unless the flow adheres exactly to the simple so- 
lution found in analytic theory, without any perturbations, 
then the required velocity nudge has to be communicated 
through shocks. 

Indeed it is this vertical velocity nudge, required twice 
per orbit, which gives rise to the compressive, pulsatory flow 
of the form Vz oc z with frequency 2Q. Again for this mo- 



tion, unless the flow adheres exactly to the analaytical solu- 
tion, and there are no perturbations present within the disc, 
then once the pulsatory motions approach sonic velocities, 
the pulsations must result in shock formation and enhanced 
dissipation. 



5 CONCLUSIONS 

Using numerical simulations we have explored the propaga- 
tion and damping of warped accretion discs in the regime 
(a ^ Qc) when propagation occurs diffusively. In this regime 
we find that warp propagation is described reasonably well 
by the simple equation given by Pringle (1992), which in- 
volves just two viscosity parameters a and 02- This is true, 
and gives the same 02 — a relationship, whether or not the 
SPH 'viscous switch' is turned on. As expected from linear 
analysis, and from the simple physical arguments presented 
in Section 4, we find that as a decreases, 02 increases and at 
large values of a our simulations agree with the predictions 
of previous analytic work. This is the first time that this has 
been demonstrated numerically. We also find, however, that, 
if a decreases below some value Osat ~ 0.16, enhanced dissi- 
pation occurs which acts to damp the warp-driven, internal, 
resonant motions, and so prevents a further rise in 02. We 
find a maximum value of Q2 — Omax ~ 3 - 4. We specu- 
late that this maximum value comes about because at low 
values of a the induced horizontal shearing motions become 
close to sonic, and are therefore subject to instabilities and 
enhanced dissipation. Further numerical work is required to 
test these ideas. 

In line with the expectations of linear theory we also 
find that the warp induces precession which causes the disc 
to twist. As predicted, the precession timescales are typi- 
cally smaller than the timescales for warp propagation and 
damping. However, the exact values of the precession rates, 
and even the direction of precession, are not in agreement 
with the predictions of linear theory. This presumably re- 
fiects at some level the degree to which SPH viscosity fails 
to provide an accurate model of isotropic Navier-Stokes vis- 
cosity, which is assumed in the linear calculations. Indeed it 
is for the simulations with the 'viscous switch' turned on (so 
that dissipation only occurs in regions of converging flow), 
for which the SPH viscosity can be expected to be least like 
Navier-Stokes, that the direction of precession is the oppo- 
site to that predicted by linear theory. Since the viscosity in 
real accretion discs is unlikely to be Navier-Stokes, this im- 
plies that theoretical predictions of warp-induced precession 
rates need to be treated with caution. 

The numerical method we have used is SPH. This has 
the advantage of being a Lagrangian code, and so can deal 
straightforwardly with thin warped discs, without the need 
of careful treatment of strong density contrasts across neigh- 
bouring grid cells. It has, however, the disadvantage that 
the fluid viscosity (a) depends on particle density, which 
therefore implies that varying viscosity involves varying nu- 
merical resolution. As can be seen from Table 1, even using 
one or two million particles, for the parameter ranges we 
are able to investigate the number of smoothing lengths per 
disc scaleheight is never very large, varying from around 10 
for the thickest disc down to 1.7 for the thinnest disc. It is 
further worth noting that although we do find evidence of 
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Figure 14. Velocity along the local direction of the angular momentum vector, in units of the local sound speed, for simulation S9, at 
t = 185, -R = 5 and = (upper left), = O.Stt (upper right), <j> = tt/2 (middle left), cf> = O.Gtt (middle right), 4> = O.Ttt (lower left) 
and at = TT (lower right). Here, the coordinate system is such that the ^-axis lies along the local direction of the angular momentum 
1. Between (f) = and </< = 7r/2 the disc undergoes supersonic compression {dvz/dz < 0), with a Mach number at z = ±H = ±0.0665 
approaching 3. Between = 0.57r and if) = 0.67r these velocities are abruptly halted. At later phases the disc undergoes supersonic 
rarefaction {dvz/dz > 0), until it reaches its original structure at </> = 7r. The sequence is repeated between <p = it and (j> = 2it because 
the orbiting disc gas has to receive a vertical velocity nudge twice per orbit in order to prevent the orbits of gaseous disc elements from 
crossing. 



enhanced dissipation as a decreases, i.e. as the amphtude 
of the internal resonant shearing motions increases, it is not 
clear that at this resolution the fluid dynamical effects which 
give rise to such enhanced dissipation (e.g. fluid instabilities, 
shocks, etc.) are being captured correctly. Thus the values 
we find for Qmax need to be treated with appropriate cau- 
tion. In addition we find that the small precessional effects 
induced by the warp depend sensitively on the exact im- 
plementation of the SPH viscosity. It would be instructive 
to compare the results we obtain here with results obtained 
using a grid based numerical technique which would be bet- 
ter able to deal with shock-capturing, and thus better able 
to identify the source of the enhanced dissipation we find at 
large amplitudes of the horizontal shear. It is also important 
to incorporate magnetic effects, since it is these which are 



thought to be the basic cause of disc 'viscosity', and it seems 
likely that these may not be well modeled by the assumption 
of an isotropic Navier-Stokes viscosity. 
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